clc;
warning('off');
Nr=10;
M1=7;
M2=9;
K=2;
R=0.5;
alpha=0.03;
e=zeros(1,K);
x_axis=1:1:5;
y_axis1=zeros(1,length(x_axis));
y_axis2=zeros(1,length(x_axis));
y_axis_random_phase1=zeros(1,length(x_axis));
y_axis_random_phase2=zeros(1,length(x_axis));
% y_axis_SDR=zeros(1,length(x_axis));
num=60;
a=0;
for i=1:num
    sign=1;
    y_axis_random_phase_temp1=zeros(1,length(x_axis));
    y_axis_random_phase_temp2=zeros(1,length(x_axis));
    y_axis_temp1=zeros(1,length(x_axis));
    y_axis_temp2=zeros(1,length(x_axis));
    disp(i);
    hr=(randn(1,Nr)+1i*randn(1,Nr))/sqrt(2);
    G1=(randn(Nr,M1)+1i*randn(Nr,M1))/sqrt(2);
    G2=(randn(Nr,M2)+1i*randn(Nr,M2))/sqrt(2);
    for Ne=x_axis
        Dr=(randn(Ne,Nr,K)+1i*randn(Ne,Nr,K))/sqrt(2);
        disp(Ne);
        for k=1:K
            e(1,k)=alpha*sqrt(norm(Dr(:,:,k),2)^2);
        end
       try
        [ w,F,power1 ] = Alternate_optimization( hr,G1,Dr,e,R,M1,Nr,Ne,K );
        [ w,F,power2 ] = Alternate_optimization( hr,G2,Dr,e,R,M2,Nr,Ne,K );
        [ power_random1 ] = random_phase( hr,G1,Dr,e,R,M1,Nr,Ne,K );
        [ power_random2 ] = random_phase( hr,G2,Dr,e,R,M2,Nr,Ne,K );
        if isnan(power_random1)
            a=a+1;
            sign=0;
            disp('false')
            break;
        end
        if isnan(power_random2)
            a=a+1;
            sign=0;
            disp('false')
            break;
        end
        catch
            a=a+1;
            sign=0;
            disp('false')
            break;
       end
        y_axis_temp1(1,int8(Ne))=power1;
        y_axis_temp2(1,int8(Ne))=power2;
        y_axis_random_phase_temp1(1,int8(Ne))=power_random1;
        y_axis_random_phase_temp2(1,int8(Ne))=power_random2;
    end
   if sign
       y_axis1=y_axis1+y_axis_temp1;
       y_axis2=y_axis2+y_axis_temp2;
       y_axis_random_phase1=y_axis_random_phase1+y_axis_random_phase_temp1;
       y_axis_random_phase2=y_axis_random_phase2+y_axis_random_phase_temp2;
   end
end
y_axis1=y_axis1/(num-a)
y_axis2=y_axis2/(num-a)
y_axis_random_phase1=y_axis_random_phase1/(num-a)
y_axis_random_phase2=y_axis_random_phase2/(num-a)
y_axis1=10*log10(y_axis1*1000)
y_axis2=10*log10(y_axis2*1000)
y_axis_random_phase1=10*log10(y_axis_random_phase1*1000)
y_axis_random_phase2=10*log10(y_axis_random_phase2*1000)
plot(x_axis,y_axis1,'b-h','LineWidth',1.5,'MarkerSize',6);
hold on;
plot(x_axis,y_axis2,'b-p','LineWidth',1.5,'MarkerSize',6);
hold on;
plot(x_axis,y_axis_random_phase1,'r-o','LineWidth',1.5,'MarkerSize',6);
hold on;
plot(x_axis,y_axis_random_phase2,'r->','LineWidth',1.5,'MarkerSize',6);
grid on;
xlabel('Number of antennas for eavesdroppers'); 
ylabel('Worst-case Transmission power (dbm)');
legend('Alternate_optimization M=7','Alternate_optimization M=9','Fix Phase M=7','Fix Phase M=9')
a